Single-ensemble nonequilibrium path-sampling estimates of free energy differences 
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We introduce a straightforward, single-ensemble, path sampling approach to calculate free en- 
ergy differences based on Jarzynski's relation. For a two-dimensional "toy" test system, the new 
(minimally optimized) method performs roughly one hundred times faster than either optimized 
"traditional" Jarzynski calculations or conventional thermodynamic integration. The simplicity of 
the underlying formalism suggests the approach will find broad applicability in molecular systems. 
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The estimation of free energy differences AF in molec- 
ular systemsi*2iMi5i£A2i£ is important for a wide variety 
of applications including virtual screening for drug de- 
sign, determination of the solubility of small molecules, 
and binding affinities of ligands to proteinsiS*ii. Jarzyn- 
ski recently introduced a general non-equilibrium ap- 
proach to computing AF 4 * 1 ^, but the technique never has 
been shown superior to more traditional equilibrium cal- 
culations (e.g. Refsi&£). Here, we introduce a potential 
route for dramatically faster non-equilibrium AF calcu- 
lations. 

Many previous workers have attempted to improve 
non-equilibrium AF estimates. Hummer studied the op- 
timization of non-equilibrium simulation^, and Jarzyn- 
ski introduced "targeted free energy perturbation" 
to improve configurational sampling^. Improvement 
of configurational sampling in AF calculations has 
also been the focus of studies by McCammon and 
collaborators^, Karplus and collaborators 1 ^ and van 
Gunsteren and collaborators^. Schulten and collabora- 
tors used Jarzynski's approach for steered molecular dy- 
namics simulations^. Ytreberg and Zuckerman 1 ^, and 
Zuckerman and Wool6*i& have developed methods for 
more efficient use of non-equilibrium data for AF calcu- 
lation. 

In an important advance of direct relevance to the 
present report, Sun suggested the use of a path sampling 
approach to evaluate AF via Jarzynski's relation, with 
a formalism that essentially entails thermodynamic in- 
tegration in (inverse) temperature spaced. Sun reported 
impressive efficiency gains. However, multiple path sam- 
pling ensembles were required even for simple systems. 

The approach outlined below builds on several sources. 
Jarzynski defined the non-equilibrium approach 4 , and 
Pratt introduced the seminal concept of sampling dy- 
namic paths with equilibrium toolset Chandler and col- 
laborators supplied Monte Carlo path sampling moves 
for effective implementation of the Pratt approachi2*i2*21i, 
and Sun suggested that path sampling ensembles could 
be used to evaluate the Jarzynski relation^. Finally, 
Zuckerman and Woolf employed a direct formalism for 
path-based estimates of arbitrary quantities, which is key 
to our single-ensemble protocol^!. 

In outline, this report first sketches Jarzynski's relation 
and shows how it can be re-written using importance 



sampling of paths. The path sampling procedure used 
in our method is then described. Finally, we present our 
results and a discussion. 

Following the usual formalism to define the AF cal- 
culation, we consider two systems or distinct states that 
are defined by Hamiltonians Hq(x) and H\{x), where x 
is a set of configurational coordinates. By introducing a 
parameter A, a hybrid Hamiltonian can be constructed, 
e.g., H(X;x) = Hq(x) + X[Hi(x) — H (x)]. Jarzynski 
showed that arbitrarily rapid, non-equilibrium switches 
from A = to A = 1 can be used to calculate the equilib- 
rium free energy difference AF = AFx=o— »i- To this end, 
one considers switching trajectories that combine incre- 
ments in A with "traditional" dynamics (such as Monte 
Carlo or Langevin dynamics) in x-space at fixed A values. 
Thus, a trajectory with n A-steps is given by 

Z„ = |(A = 0,x ), (Ai,fo), (Ai,fi), (A 2 ,xi), 

(A 2 ,x 2 ), (Xn-l,X n -i), (A„ = l,x„_i)|, (1) 

where it should be noted that increments (steps) from Aj 
to Ai+i are performed at a fixed conformation Xi, and 
the initial xq is drawn from the Hq distribution. For 
simplicity we have assumed only a single dynamics step at 
each fixed Xi, from to Xi, is performed, but multiple 
steps can be performed within the Jarzynski formalism. 

Finally, the work performed on the system during a 
switching trajectory is 



n-l 

W r (Z„) = ^[i?(A i+1 ;f i )-lf(A i 



(2) 



and transcribing the Jarzynski relation* into path lan- 
guage, the free energy difference can be written as^ 



-/3AF _ 



J dZ n Q(Z n ) e -^ z "> / JdZ n Q(Z„), 



(3) 



where [3 = 1/ksT, dZi n denotes integration over all pos- 
sible trajectories, and Q(Z n ) is proportional to the prob- 
ability of occurrence of trajectory Z„. Q(Z n ) depends 
on the dynamics employed and will be specified below 
for the overdamped Langevin case. 
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In "standard" non-equilibrium simulation, the integral 
P|) need never be considered since trajectories — and the 
associated work values — are automatically generated 
with the proper frequency (i.e., proportional to Q(Z n )). 
In this case, the Jarzynski relation provides an estimate 
for AF for a set of work values {W\, W2, Wn} given 
byi 



AF = AF Jarz = -~ln 
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(4) 



where the "=" denotes a computational estimate. Since 
the relationships in Eqs. Q and (0J are valid for an arbi- 
trary number n of A-steps, switches may be performed 
very rapidly. The apparent advantage of these "fast- 
growth" (small n) calculations is that very little com- 
putational time is spent generating trajectories, and thus 
AFj arz can be generated with very little CPU time. How- 
ever, in practice, unless there is sufficient overlap between 
the states described by H (x) and Hi(x), AFj arz will be 
biased, often by many fcg T 5 ' 6 ' 13 . This bias is due to the 
nonlinear nature of Eq. (0} where the smallest, and thus 
rarest, work values dominate the average. Additionally, 
CPU time must be invested in generating the equilibrium 
distribution for Hq. 

This study uses importance sampling of switching tra- 
jectories to sample dominant but rare work values more 
frequently, without the need to sample the Hq equilib- 
rium distribution. We combine the sampling strategy 
of Sun£ with the simple formalism used by Zuckerman 
and WoolfSi, as we consider an alternative distribution 
of switching trajectories D{Z n ). Then, with no loss of 
generality, Eq. @ can be written as 



-/3AF 



/ dZ n D(Z n ) 


Q(Z n )/£>(Z„) 




J dZ n D(Z n ) 


Q(Z n )/D(Z n ) 





E D Q(z»)e-^( z ")/j(z») 

£ D Q(Z„)AD(Z„) 



(5) 



where the only condition is that D(Z n ) ^ anywhere. 
The shorthand indicates a sum over trajectories gen- 
erated according to D(Z n ). 

Since the fundamental idea behind the importance 
sampling in Eq. (JSJ is to generate trajectories — and 
hence work values — according to D(Z n ), the choice D 
is critical. We choose D(Z n ) to favor trajectories with 
important work values, namely, 



D(Z n ) = Q(Z n )e-^ w ^. 



(6) 



As will be seen below in Eq. J7J), this choice appears to 
balance convergence difficulties between the numerator 
and denominator of Eq. JSJ). We note that Sun also em- 
ployed the distribution (JSJ) as one among several used for 
an indirect calculation of AF— . While it is not obvious 
that the choice is optimal in general, other forms for 



D(Z n ) have been tested by the authors and provided no 
improvement. By comparison with 0, the (3/2 in © 
embodies double the temperature. 

Combining Eqs. (JSJ and ©, the free energy estimate 
for our single-ensemble path sampling (SEPS) method is 
given by the new relation 



AF = AF sops = --ln 



■ hpw 



(7) 



We now specify Q(Z n ) from Eq. J7J) which is required 
for the path sampling performed below. We assume over- 
damped Langevin (Brownian) dynamics is used at fixed A 
values. Single-step distributions for ASi = Xi — Xi— 1 are 
thus Gaussian, with a variance given by a 2 = 2At/m-y(3, 
where m is the mass of the particle and 7 is the friction 
coefficient of the medium (e.g. RefiSi). Combining the 
Brownian distributions with that for A = leads to the 
full trajectory weight 



a ~\Ax i -Axf Bt \ 2 /2o 

Q(Z n )=e-^^n 

i=l 



(27R7 2 ) d / 2 



(8) 



where Axf ct — —V x H(Ai;xi^i)At/mj is proportional 
to the force and time step, and d is the dimensionality 
of the conformational space x. Note that if deterministic 
dynamics (e.g., Verlet) is used then Q(Z n ) = e~P H °( x °). 

To calculate the free energy estimate AF scps in Eq. 
(|7|). switching trajectories must be generated according 
to D(Z n ). This is readily accomplished using the path 
sampling approach proposed by Prattii, where entire tra- 
jectories (paths) are generated and then accepted or re- 
jected based upon a suitable Monte Carlo criteria. Trial 
moves in path space are generated following Chandler 
and coworkersi&iii. 

Putting the pieces together, we estimate AF by sam- 
pling trajectories according to D(Z n ) in Eq. © using 
the following steps (c.f. Ref. 3 ): (i) Generate an arbi- 
trary initial reference trajectory by switching the sys- 
tem from A = — > 1. Calculate the work W done 
on the system during the switch, (ii) Pick a random 
A value along the reference trajectory and make a ran- 
dom phase-space displacement. For Brownian dynamics 
this corresponds to a random shift in position. Gen- 
erate a trial trajectory by "shooting" forward (incre- 
ment A) and backward (decrement A). Calculate the 
trial work done on the system W'. (iii) Accept this 
new trajectory according to the Metropolis criterion: 
mm[l,Q'e-il 3W /Qe-§P W ], with Q from Eq. ©■ (iv) If 
accepted, the trial trajectory becomes the current refer- 
ence trajectory. If rejected, the current reference trajec- 
tory remains unchanged. Whether accepted or rejected, 
the current reference trajectory is then used in Eq. J7J). 
Repeat from step (ii). 

It should be noted that to obtain good sampling, as 
in any Monte Carlo simulation, equilibrium must be at- 
tained before averages are calculated. Using the path 
sampling procedure above, we accomplish this by check- 
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FIG. 1: Contour plots of the test system Ho(x,y) and 
Hi(x,y) given by Eq. @ where each contour represents an 
energy change of 4.0 fcsT. This problem is expected to be dif- 
ficult due to the large barrier and the asymmetric double- well 
in Hi. 



ing the running average work every 20 accepted trajec- 
tories for convergence within 0.01 fcgT. 

As a test problem, we consider a two-dimensional sys- 
tem (Fig. ^) switched from a single well to a double well: 

H (x,y) = (x + 2) 2 + y 2 , 

H x {x,y) = !{((* - I) 2 - y 2 ) 2 + 10(x 2 - 5) 2 + 

(x + y)* + (x-y)*}. (9) 

Figurendearly demonstrates why estimating AF for this 
system is expected to be difficult: the significant barrier 
in Hi will prevent sufficient configurational sampling of 
the dominant minimum at -Hi (2,0) for short trajecto- 
ries. Thus ordinary fast-growth Jarzynski estimates will 
substantially overestimate AF. Similarly, equilibrium 
approaches like thermodynamic integration (TI) will re- 
quire long simulation times to surmount the barrier. 

For this system, the free energy difference was esti- 
mated using the Jarzynski method given by AFj arz in 
Eq. (QJ, by the SEPS method given by AF scps in Eq. 
0, as well as by conventional TI. Trajectories for all es- 
timates were generated using Brownian dynamics with 
parameters j3 = 7 = m = 1, and At = 0.001. 

For AFj arz , to generate uncorrelated initial configura- 
tions xq, the system was run at A = for N cq steps be- 
tween switching trajectories. For Hq given by @, it was 
determined that N cq = 10,000, and that smaller values 
of iVeq introduce bias in AFj arz . Given N eq , moreover, 
we optimized AFj arz by varying the number n of X-steps 
in Eqs. (JTJ and ©• We found that n = 100,000 was 
most efficient. 

Trajectories for AF sops were generated as described 
above. Specifically, perturbations to the selected state 
Xi of the reference trajectory (step (ii) above) were cho- 
sen from a Gaussian distribution of width 50.0a, giving 
an acceptance ratio of 1 — 2%. Smaller perturbations 
were also highly successful. The SEPS procedure is not 
optimized in the sense that only a simple type of trial 
move (termed "shooting"—) was employed, and we used 
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FIG. 2: Comparison between free energy estimates from the 
Jarzynski method, conventional Thermodynamic integration 
(TI), and our single-ensemble path sampling (SEPS) method. 
The circles show the results of the SEPS method using 10 
A-steps. The results of the Jarzynski method for a very short 
trajectory (10 A-steps, squares) and the most efficient trajec- 
tory length (100,000 A-steps, triangles) are also shown. TI 
estimates based on 10 A increments are shown as diamonds. 
The exact answer of AF = 6.55 fcsT is shown as a solid hor- 
izontal line. Each data point represents the mean estimate, 
with standard deviations given by the error bars, based on 
100 independent estimates of AF for each method. 



strict path-equilibration criteria. Optimization methods 
are currently under investigation by the authors. 

For comparison to an equilibrium approach, "text- 
book" thermodynamic integration (TI) simulations 22 
were performed with identical Brownian parameters and 
10 A-steps, with 25% of data discarded for equilibration. 
This well-known approach is described in many sources 
(e.g., Ref^) and is not detailed here. Since no optimiza- 
tion was performed, wc refer to this method as "conven- 
tional TI." 

To compare the efficiency of the SEPS approach with 
other methods, in Fig. [21 we plot AF estimates for the 
SEPS, Jarzynski, and TI methods as a function of the to- 
tal CPU time needed generate the estimates. The circles 
show the results of the SEPS method using 10 A-steps. 
Also shown are the results of the Jarzynski method us- 
ing a very short trajectory (10 A-steps, squares), and the 
most efficient trajectory length (100,000 A-steps, trian- 
gles). TI estimates based on 10 A increments are shown 
as diamonds. The solid horizontal line gives the exact 
answer AF = 6.55 fc^T. The plot was generated by cal- 
culating the mean (data points) and standard deviations 
(error bars) from 100 independent estimates of AF for 
each method. The CPU time spent equilibrating is in- 
cluded in the total CPU time for all methods. 

As expected, Fig. [21 shows that for fast-growth work 
values (10 A-steps, squares), the Jarzynski method incor- 
rectly estimates the free energy difference as AFj arz w 
13ksT. As the number of A-steps increases, the stan- 
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dard Jarzynski trajectories begin to "see" the minimum 
at £Zi(2,0) and the correct AF is obtained. Since the 
highest efficiency for the standard Jarzynski method was 
obtained using f 00, 000 A-steps (triangles), we consider 
this curve to be the optimized Jarzynski method for the 
test system. The unoptimized, conventional TI calcu- 
lations are of comparable efficiency to the traditional 
Jarzynski estimates. 

The SEPS method, by contrast, correctly estimates 
the free energy quickly and accurately, even for very 
short trajectories (10 A-steps). One can quantitatively 
compare estimates from each method by noting that the 
estimate for the SEPS method AF seps (f w 1500 s) is 
slightly more accurate than Ai 7 j arz (t w 150,000 s) and 
AF T i(t « 500, 000 s), implying a roughly 100-fold speed- 
up of SEPS over the other methods. 

Compared to "standard" Jarzynski calculation, the 
SEPS approach has several advantages: (a) important, 
rare trajectories with small work values are favored; (b) 
no CPU time is spent acquiring an equilibrium ensemble 
at A = 0; and (c) path-sampling moves that are capable 
of surmounting barriers may be used. In other words, 
the SEPS approach focuses CPU time on the important 
regions of (A; x) space — this also contrasts with TI and 



other equilibrium approaches which attempt to sample 
the full space. 

To summarize, we have described a rapid and straight- 
forward new method for estimating free energy differ- 
ences AF, using a single-ensemble path sampling (SEPS) 
approach. We also have carefully quantified the numer- 
ical efficiency of the approach. Without extensive op- 
timization, the SEPS method generates AF estimates 
over 100 x more efficiently than "standard" Jarzynski and 
conventional thermodynamic integration calculations, for 
the two-dimensional test system considered here. Our 
approach relies on an extremely simple importance sam- 
pling formalism, and therefore appears to be readily ex- 
tendable to molecular systems. This extension — which 
will require addressing issues of memory for trajectory 
storage — is currently underway. We will also optimize 
the SEPS approach via alternative importance sampling 
distributions, and path-sampling trial moves. 
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